OpenCyto Analysis of HVTN080

The HVTN080 study flow cytometry data has been preprocessed (loaded in R from the raw FCS files and the FlowJo workspaces). The gating set is available here, and has been gated using OpenCyto.

We have gated the data using this OpenCyto template.

The manually gated data is available here, and contains the manual gates from the FlowJo workspaces.

The raw FCS files are available from flowrepository.org/id/FR-FCM-ZZ7U.

The data in the gating sets has already been compensated and transformed using the compensation and transformation information stored in the FlowJo workspaces.

The FlowJo workspaces for this data will be posted soon.

The code to perform the gating and extract the features is here.

Visualization of the gating hierarchy

We load the gated data and visualize the results.

suppressPackageStartupMessages(
  {require(flowWorkspace)
  require(flowViz)
  require(knitr)
  require(data.table)
  require(MIMOSA)
  require(GFMisc)
  require(scales)
  require(lme4)
  require(contrast)
  require(multcomp)
  require(COMPASS)
})
opts_chunk$set(list(message=FALSE,warning=FALSE))
WORKDIR<-"/Users/gfinak/Dropbox/GoTeam/Projects/Paper-OpenCytoPipeline/HVTN080/reproducible/"
AUTO<-"./gs_auto_clean"
MANUAL<-"./gs_manual_clean"
setwd(WORKDIR)
auto<-load_gs(AUTO)

## loading R object...
## loading tree object...
## Done

manual<-load_gs(MANUAL)

## loading R object...
## loading tree object...
## Done

Here is a view of the gating hierarchy for the manually gated data.

plot(manual[[1]])

plot of chunk unnamed-chunk-3

And here is the same for the OpenCyto gated data.

plot(auto[[1]])

plot of chunk unnamed-chunk-5

Visualize gating layouts for manual and automated gating

We can view dotplots of each gate from selected samples for each of the OpenCyto and manually gated data.

First, the OpenCyto gated data:

plotGate(auto[[1]],scales=list(cex=c(1.2,1.2)),par.strip.text=list(cex=1.4),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=1)),path=2)

plot of chunk unnamed-chunk-6

Then the manually gated data:

plotGate(manual[[1]],scales=list(cex=c(1.2,1.2)),par.strip.text=list(cex=1.4),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=1)),path=2)

Comparison of the perforin gate

We want to look at the variablity of the perforin channel in the OpenCyto analysed data.

#extract manually gated perforin for CD8+ for a single sample
#p1<-plotGate(manual["505325.fcs"],"8+/Perforin+",scales=list(cex=c(1.2,1.2)),par.strip.text=list(cex=1.1),xlab=list(fontsize=18),ylab=list(fontsize=18),par.settings=list(gate.text=list(cex=2)),isPath=FALSE,main="Manual Gate",arrange=FALSE)

#Extract the same gate from the OpenCyto gated data
set.seed(100)
s<-sample(1:length(auto),6)
p2<-plotGate(auto[s],"8+/Perforin+",scales=list(cex=c(1.1,1.1)),par.strip.text=list(cex=1.1),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=1)),isPath=FALSE,default.y="<Alexa 680-A>",main="OpenCyto Gate",arrange=FALSE)

#generate a nice layout in a single panel
pdf(file="../../Submission/Figure4.pdf",width=8,height=5)
grid.arrange(p2,main="Perforin Expression")
dev.off()

## pdf 
##   2

tiff(file="../../Submission/Figure4.tiff",width=8*150,height=5*150,res=150)
grid.arrange(p2,main="Perforin Expression")
dev.off()

## pdf 
##   2

grid.arrange(p2,main="Perforin Expression")

plot of chunk unnamed-chunk-8

Side-by-side comparison of OpenCyto and manual gating for other gates

We can generate any layout we want, even viewing plots for the same gate and sample side-by-side from the manual and openCyto analysis.

We can control the display projections of individual gates.

# Comparison of manual vs Autogating for ICS
proj=list("S"=c(x="FSC-H",y="FSC-A"),"L"=c(x="FSC-A",y="SSC-A"),"4+"=c(x="CD8",y="CD4"),"3+"=c(y="CD3",x="CD8"),"8+/57+"=c(x="CD57",y="Granzyme B"),"8+/Granzyme B+"=c(x="Perforin",y="Granzyme B"),"4+/IL2+"=c(x="IL2",y="IFNg"),"4+/IFNg+"=c(x="IL2",y="IFNg"),"8+/IL2+"=c(x="IL2",y="IFNg"),"8+/IFNg+"=c(x="IL2",y="IFNg"))
man<-plotGate(manual[["505325.fcs"]],c("S","Lv","L","3+","4+","8+","4+/IFNg+","4+/IL2+","8+/IFNg+","8+/IL2+","8+/Granzyme B+","8+/57+"),projections=proj,scales=list(cex=c(0.7,0.7)),par.strip.text=list(cex=0.7),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=0)),path=1,arrange=FALSE,xbin=64,merge=TRUE)
aut<-plotGate(auto[["505325.fcs"]],c("S","Lv","L","3+","4+","8+","4+/IFNg+","4+/IL2+","8+/IFNg+","8+/IL2+","8+/Granzyme B+","8+/57+"),projections=proj,scales=list(cex=c(0.7,0.7)),par.strip.text=list(cex=0.7),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=0)),path=1,arrange=FALSE,xbin=64,merge=TRUE)
gates<-list()
for(i in 1:9){
  gates[[i]]<-c(man[[i]],aut[[i]])
}
grid.arrange(do.call(arrangeGrob,gates),ncol=1)

plot of chunk unnamed-chunk-9

pdf(file="../../Submission/Figure2.pdf",width=8,height=6)
grid.arrange(do.call(arrangeGrob,gates),ncol=1)
dev.off()

## pdf 
##   2

tiff(file="../../Submission/Figure2.tiff",width=8*150,height=6*150,res=150)
grid.arrange(do.call(arrangeGrob,gates),ncol=1)
dev.off()

## pdf 
##   2

Evaluation of vaccine-specific effects in Perforin-positive subsets of cells

Since OpenCyto was able to accurately gate perforin, we’d like to see if we can detect any vaccine-specific perforin cell subsets.

We extract cell counts from the CD4 and CD8 T-cell subsets that express perforin in combination with other cytokines, using COMPASS.

First, we extract the counts from the automated gates using COMPASS. We focus on cells expressing any cytokine together with perforin, and test whether the proportion of pre vs post-vaccine cells differs by treatment group (Wilcoxon test) within the CD4 and CD8 subsets.

#extract the population statistics
auto_stats<-getPopStats(auto,statistic="count")
if(!file.exists("./CD8_CompassContainer.rds")){
cd8_auto<-COMPASSContainerFromGatingSet(gs=auto,node="8\\+$",individual_id="PTID",sample_id="name",markers=c("IFNg","TNFa","IL2","Granzyme B","Perforin"))
saveRDS(cd8_auto,file="/Users/gfinak/Dropbox/GoTeam/Projects/Paper-OpenCytoPipeline/HVTN080/CD8_CompassContainer.rds")
}else{
  cd8_auto<-readRDS(file="./CD8_CompassContainer.rds")
}
if(!file.exists("./CD4_CompassContainer.rds")){
cd4_auto<-COMPASSContainerFromGatingSet(gs=auto[pData(auto)$Stim%in%c("POL-1-PTEG","negctrl 1","negctrl2 ")],node="4\\+$",individual_id="PTID",sample_id="name",markers=c("IFNg","TNFa","IL2","Granzyme B","Perforin"))
saveRDS(cd4_auto,file="/Users/gfinak/Dropbox/GoTeam/Projects/Paper-OpenCytoPipeline/HVTN080/CD4_CompassContainer.rds")
}else{
  cd4_auto<-readRDS("./CD4_CompassContainer.rds")
}

#Plots the proprotions for a given combination for CD8
.plotProps<-function(cd8,title,combination){
  #Compute cell counts for the given combination
  C<-COMPASS:::CellCounts(cd8,combinations=combination)
  tofit<-cbind(cd8$meta,cbind(cd8$counts-C,C))
  colnames(tofit)[c(25,26)]<-c("NSUB","CYTNUM")
  tofit$Stim<-factor(tofit$Stim)
  #Aggregate negative controls
  levels(tofit$Stim)<-c("ENV-1-PTEG","GAG-1-PTEG","negctrl","negctrl","POL-1-PTEG")
  
  
  #relevel the visit variable
  levels(tofit$VISITNO)<-c("Pre-vaccine","Post-vaccine")
  tofit$rx_code<-factor(tofit$rx_code)
  
  #relevel the treatment variable
  levels(tofit$rx_code)<-c("Placebo","Pennvax B", "Pennvax B + IL12 DNA")
  
  tofit2<-ddply(tofit,.(Stim,VISITNO,PTID),function(x)with(x,data.frame(rx_code=unique(rx_code),Stim=unique(Stim),CYTNUM=sum(CYTNUM),NSUB=sum(NSUB))))
  p1<-ggplot(subset(tofit2,Stim%in%"POL-1-PTEG"))+geom_boxplot(aes(x=rx_code,y=CYTNUM/(NSUB+CYTNUM)))+facet_grid(Stim~VISITNO)+theme_bw()+scale_y_continuous("Empirical Proportion of Vaccine\nResponsive Cells")+scale_x_discrete("Treatment Group")+coord_trans(ytrans=log1p_trans())
  E<-ConstructMIMOSAExpressionSet(tofit2,reference=Stim=="negctrl",measure.columns=c("NSUB","CYTNUM"),other.annotations=c("rx_code","VISITNO","Stim","PTID"),default.cast.formula=component~PTID+Stim+VISITNO+rx_code,.variables=.(PTID,rx_code,VISITNO))
  fitted<-MIMOSA(NSUB+CYTNUM~rx_code+PTID+RefTreat+Stim+VISITNO,E,ref=RefTreat%in%"Reference"&Stim%in%"POL-1-PTEG",subset=RefTreat%in%"Treatment"&Stim%in%"POL-1-PTEG",method="mcmc",pXi=1,EXPRATE=1e-6,tune=500)

  #linear model contrast pr vs post vaccine within treatment groups  
  lmdata<-cbind(with(data.frame(countsTable(fitted,proportion=TRUE)),data.frame(proportion=CYTNUM-CYTNUM_REF)),"Pr.Response"=getZ(fitted)[,2],pData(fitted))
  KK<-GFMisc:::optimAsinhCofactor(subset(lmdata,rx_code%in%"Placebo")$proportion)
  lmdata$prop_trans<-asinh(lmdata$proportion*KK)
  lmfit<-lm(prop_trans~rx_code*VISITNO,data=lmdata)
  contr<-contrast(lmfit
                  ,a=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Post-vaccine"),
                  b=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Pre-vaccine"))
  fit1<-(lmer(prop_trans~rx_code*VISITNO+(1|PTID),data=lmdata))
  sumry<-data.frame(p.value=summary(glht(fit1,linfct=contr$X,alternative="greater"))$test$pvalues,rx_code=c("Pennvax B","Pennvax B + IL12 DNA")) 
  lmdata<-data.table(lmdata)
  setkey(lmdata,PTID,VISITNO)
  #extract the data and plot
  lmdata_plot<-lmdata[,list(difference=proportion[2]-proportion[1]),list(rx_code,PTID)]
  p<-ggplot(subset(lmdata_plot,!rx_code%in%"Placebo"),aes(x=rx_code))+geom_boxplot(aes(y=difference))+theme_bw()+geom_text(aes(y=3e-4,label=paste0("p = ",signif(p.value,3))),data=sumry)+scale_y_continuous("Post-Vaccine - Pre-Vaccine")+scale_x_discrete("Treatment Group")
  
  return(p)
  }

#Plots the proportion for a given combination for CD4
.plotPropsCD4<-function(cd8,title,combination){
  C<-COMPASS:::CellCounts(cd8,combinations=combination)
  tofit<-cbind(cd8$meta,cbind(cd8$counts-C,C))
  colnames(tofit)[c(25,26)]<-c("NSUB","CYTNUM") 
  tofit$Stim<-factor(tofit$Stim)
  
  #Aggregate negative controls
  levels(tofit$Stim)<-c("negctrl","POL-1-PTEG")
  
  
  #relevel the pre-post vaccine variable
  levels(tofit$VISITNO)<-c("Pre-vaccine","Post-vaccine")
  tofit$rx_code<-factor(tofit$rx_code)
  
  #relevel the treatment groups
  levels(tofit$rx_code)<-c("Placebo","Pennvax B", "Pennvax B + IL12 DNA")
  
  #some data reshaping
  tofit2<-ddply(tofit,.(Stim,VISITNO,PTID),function(x)with(x,data.frame(rx_code=unique(rx_code),Stim=unique(Stim),CYTNUM=sum(CYTNUM),NSUB=sum(NSUB))))
  E<-ConstructMIMOSAExpressionSet(tofit2,reference=Stim=="negctrl",measure.columns=c("NSUB","CYTNUM"),other.annotations=c("rx_code","VISITNO","Stim","PTID"),default.cast.formula=component~PTID+Stim+VISITNO+rx_code,.variables=.(PTID,rx_code,VISITNO))
  # We don't need to fit a MIMOSA model here, I just want to extract the cell counts with the metadata.
  fitted<-MIMOSA(NSUB+CYTNUM~rx_code+PTID+RefTreat+Stim+VISITNO,E,ref=RefTreat%in%"Reference"&Stim%in%"POL-1-PTEG",subset=RefTreat%in%"Treatment"&Stim%in%"POL-1-PTEG",method="mcmc",pXi=1,EXPRATE=1e-6,tune=500)
  
  
  #linear model testing pre vs post-vaccine within each treatment group
  lmdata<-cbind(with(data.frame(countsTable(fitted,proportion=TRUE)),data.frame(proportion=CYTNUM-CYTNUM_REF)),"Pr.Response"=getZ(fitted)[,2],pData(fitted))
  KK<-GFMisc:::optimAsinhCofactor(subset(lmdata,rx_code%in%"Placebo")$proportion)
  lmdata$prop_trans<-asinh(lmdata$proportion*KK)
  lmfit<-lm(prop_trans~rx_code*VISITNO,data=lmdata)
  contr<-contrast(lmfit
                  ,a=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Post-vaccine"),
                  b=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Pre-vaccine"))
  fit1<-(lmer(prop_trans~rx_code*VISITNO+(1|PTID),data=lmdata))
  sumry<-data.frame(p.value=summary(glht(fit1,linfct=contr$X,alternative="greater"))$test$pvalues,rx_code=c("Pennvax B","Pennvax B + IL12 DNA")) 
  lmdata<-data.table(lmdata)
  setkey(lmdata,PTID,VISITNO)
  lmdata_plot<-lmdata[,list(difference=proportion[2]-proportion[1]),list(rx_code,PTID)]
  #Pull out the results and plot
  p<-ggplot(subset(lmdata_plot,!rx_code%in%"Placebo"),aes(x=rx_code))+geom_boxplot(aes(y=difference))+theme_bw()+geom_text(aes(y=3e-4,label=paste0("p = ",signif(p.value,3))),data=sumry)+scale_y_continuous("Post-Vaccine - Pre-Vaccine")+scale_x_discrete("Treatment Group")
return(p)
}


#proportion of cells expressing any cytokine and perforin, by treatment group
cap<-capture.output(anycyt<-.plotProps(cd8=cd8_auto,title="POL-1-PTEG Stimulation of\nCD8+ T-cells for (IL2 or IFNg or TNFa) AND Perforin",combination="(IL2|IFNg|TNFa)&Perforin"))
cap<-capture.output(cd4_anycyt<-.plotPropsCD4(cd8=cd4_auto,title="POL-1-PTEG Stimulation of CD4+ T-cells\nfor (IL2 or IFNg) AND Perforin",combination="(IL2|IFNg)&Perforin"))
grid.arrange(arrangeGrob(arrangeGrob(anycyt,main="CD8 (Any Cytokine & Perforin)"),arrangeGrob(cd4_anycyt,main="CD4 (Any Cytokine & Perforin)"),nrow=1))

plot of chunk unnamed-chunk-11

SessionInfo for reproducibility

sessionInfo()

## R Under development (unstable) (2014-02-06 r64933)
## Platform: x86_64-apple-darwin13.0.0 (64-bit)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] splines   parallel  grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] sandwich_2.3-0       COMPASS_1.1.5-1      multcomp_1.3-2      
##  [4] TH.data_1.0-3        mvtnorm_0.9-9997     contrast_0.19       
##  [7] rms_4.1-3            SparseM_1.03         Hmisc_3.14-0        
## [10] Formula_1.1-1        survival_2.37-7      lme4_1.0-6          
## [13] Matrix_1.1-2         scales_0.2.3         GFMisc_1.0          
## [16] MIMOSA_0.99.5        ggplot2_0.9.3.1      Biobase_2.23.4      
## [19] BiocGenerics_0.9.3   reshape_0.8.4        plyr_1.8            
## [22] MASS_7.3-29          data.table_1.8.11    knitr_1.5.22        
## [25] flowWorkspace_3.9.65 gridExtra_0.9.1      ncdfFlow_2.9.19     
## [28] flowViz_1.27.16      lattice_0.20-27      flowCore_1.29.23    
## 
## loaded via a namespace (and not attached):
##  [1] abind_1.4-0          clue_0.3-47          cluster_1.14.4      
##  [4] coda_0.16-1          colorspace_1.2-4     corpcor_1.6.6       
##  [7] DEoptimR_1.0-1       dichromat_2.0-0      digest_0.6.4        
## [10] evaluate_0.5.1       feature_1.2.10       formatR_0.10        
## [13] graph_1.41.2         gtable_0.1.2         hexbin_1.26.3       
## [16] IDPmisc_1.1.17       KernSmooth_2.23-10   ks_1.8.13           
## [19] labeling_0.2         latticeExtra_0.6-26  MCMCpack_1.3-3      
## [22] minqa_1.2.3          modeest_2.1          munsell_0.4.2       
## [25] nlme_3.1-113         pcaPP_1.9-49         pracma_1.6.4        
## [28] proto_0.3-10         RColorBrewer_1.0-5   Rcpp_0.11.0.2       
## [31] reshape2_1.3.0.99    Rgraphviz_2.7.0      rmarkdown_0.1.5     
## [34] robustbase_0.90-2    rrcov_1.3-4          RSVGTipsDevice_1.0-4
## [37] stats4_3.1.0         stringr_0.6.2        testthat_0.8.1      
## [40] tools_3.1.0          XML_3.98-1.1         yaml_2.1.10         
## [43] zlibbioc_1.9.0       zoo_1.7-11